.
mae <- readRDS("~/BHK lab/Ravi/Ravi_Test/data/result3.rds")
#removing na in response and expression data
na_rows <- is.na(mae$response)
cleaned_response <- mae$response[!na_rows]
cleaned_expr_data <- assays(mae)[["expr"]][, !na_rows]
# 122 patient ate removing NA
design <- model.matrix(~ cleaned_response)
fit <- lmFit(cleaned_expr_data, design)
fit <- eBayes(fit)
#Top table only shows the top portion of the results
topTable(fit)
## Removing intercept from test coefficients
## logFC AveExpr t P.Value adj.P.Val
## ENSG00000211835.1 2.9729828 -7.868755 4.290500 3.620134e-05 0.9693367
## ENSG00000173231.6 1.5775015 -8.827553 3.826189 2.077302e-04 0.9693367
## ENSG00000229037.2 1.5434391 -8.867183 3.734687 2.884608e-04 0.9693367
## ENSG00000261335.1 1.0303203 -9.245205 3.693547 3.337526e-04 0.9693367
## ENSG00000251497.2 1.0942021 -9.070539 3.660129 3.754202e-04 0.9693367
## ENSG00000244053.1 1.7543381 -2.320355 3.491816 6.713542e-04 0.9693367
## ENSG00000236651.1 0.8013953 -9.597931 3.478346 7.027426e-04 0.9693367
## ENSG00000232549.1 2.0167536 -7.511261 3.469648 7.237379e-04 0.9693367
## ENSG00000232173.2 1.4258124 -8.763200 3.465522 7.339043e-04 0.9693367
## ENSG00000266707.1 2.2187534 -6.882755 3.457535 7.539641e-04 0.9693367
## B
## ENSG00000211835.1 1.85507430
## ENSG00000173231.6 0.41670113
## ENSG00000229037.2 0.14769222
## ENSG00000261335.1 0.02836268
## ENSG00000251497.2 -0.06781811
## ENSG00000244053.1 -0.54187398
## ENSG00000236651.1 -0.57905328
## ENSG00000232549.1 -0.60299910
## ENSG00000232173.2 -0.61434367
## ENSG00000266707.1 -0.63627069
datatable(topTable(fit))
## Removing intercept from test coefficients
Preparing Data for volcano Plot
Convert fit to Data frame and add column gene_name symbol
result <- topTable(fit, number=Inf) # Get all results
## Removing intercept from test coefficients
df <- as.data.frame(result) # Convert to data frame
#extract gene_name from mae ( it has gene_id column and gene_name)
genedata<- data.frame(rowData(mae@ExperimentList$expr))
#subset gene_id and gene_name from the gene_data
subset_genedata <- genedata[, c("gene_id","gene_name" )]
#mergeing process
#add a gene_id column to result
result$gene_id <- rownames(result)
#merge result and subset_genedata by gene_id
merge_result <- merge(result, subset_genedata, by= "gene_id")
datatable(merge_result)
## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html
Volcano Plot
EnhancedVolcano(merge_result,
lab = merge_result$gene_name,
x = 'logFC',
y = 'P.Value')